Nature of the spin liquid state of the Hubbard model on honeycomb lattice 
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Recent numerical work (Nature 464, 847 (2010)) indicates the existence of a spin liquid phase (SL) 
that intervenes between the antiferromagnetic and semimetallic phases of the half filled Hubbard 
model on a honeycomb lattice. To better understand the nature of this exotic phase, we study the 
quantum J 1 — J2 spin model on the honeycomb lattice, which provides an effective description of the 
Mott insulating region of the Hubbard model. Employing the variational Monte Carlo approach, we 
analyze the phase diagram of the model, finding a phase transition between antiferromagnet and an 
unusual Z2 SL state at J2/J1 ~ 0.08, which we identify as the SL phase of the Hubbard model. At 
higher J2/J1 > 0.3 we find a transition to a dimerized state with spontaneously broken rotational 
symmetry. 



Introduction. The Hubbard model describes elec- 
trons hopping on a lattice and interacting via on-site 
Coulomb interactions. 
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where s =T, J, denotes spin and m s = a is a,i s . Despite 
its conceptual simplicity, the Hubbard model exhibits a 
rich phase diagram and is believed to capture the physics 
of the high-temperature cuprate superconductors [3] (for 
a review, see Ref. [1J). At half- filling and strong repul- 
sion, (7>t, the Hubbard model is a Mott insulator, in 
which electrons are localized by strong Coulomb repul- 
sion. The Mott insulator is characterized by a charge gap 
of the order U, and, in the limit t <C U the low-energy 
dynamics of this phase is associated with the spin degree 
of freedom. The effective spin-spin interactions, which 
originate from virtual hopping processes and intertwine 
spins of neighboring electrons, become increasingly frus- 
trated as the ratio U/t is lowered and the Mott transi- 
tion is approached. In the vicinity of the Mott transition, 
the frustration enhances quantum fluctuations, which can 
prevent ordering of spins down to zero temperature, giv- 
ing rise to spin liquid ground states. The interest in spin 
liquids stems from the fact that some of them exhibit 
new types of topological order [2], and the fact that their 
properties may be linked with the physics of the doped 
Hubbard model [T]. 

Recently, the Hubbard model on the honeycomb lat- 
tice at half-filling was studied using the determinantal 
quantum Monte Carlo (DQMC) method [4 . The Hub- 
bard model on the honeycomb lattice has a crucial ad- 
vantage of being free of the sign problem, and therefore 
DQMC gives essentially exact results for correlators of 
the system. It was found that in the vicinity of the 
Mott transition, the system exhibits a disordered spin 
phase. This phase intervenes between the antiferromag- 
netic Neel state realized at higher U/t ss 4.3, and the 
semimetallic phase at U/t < 3.5. The authors of Ref. [4 
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FIG. 1: Phase diagram of the quantum J\ — J2 model, 
obtained using the variational Monte Carlo method. At 
J2/J1 ~ 0.08 we have found a phase transition between an 
antiferromagnetic state, realized in a pure Heisenberg model, 
and a gapped spin liquid state. The arrow indicates the lo- 
cation for this transition in the Hubbard model [4]. The SL 
state is best described variationally by SPS. The spin model 
at J2/J1 ^ 0.12 describes the Mott insulating phase of the 
Hubbard model. At higher J2/J1 > 0.3 the SL gives way to 
a dimerized phase. The trial wave functions of AFM states 
were of the Huse-Elser type |10| . and for SL and dimerized 
state of the RVB type |Bj. 



found that the disordered phase shows a small but fi- 
nite spin gap, and preserves translational symmetry and 
time-reversal symmetry. This suggests that the disor- 
dered state on the honeycomb lattice is a non-chiral spin 
liquid. 

In this paper, we elucidate the nature of the spin liquid 
state (SL) on the honeycomb lattice. We study the effec- 
tive J\ — J2 spin model of the large-[7 Hubbard model. 
Using Variational Monte Carlo (VMC) method, we find 
a phase diagram (FIg.l), which leads us to identify the 
exotic phase seen in the Hubbard model as the sublattice- 
pairing-state (SPS), a small-gap Z 2 SL, first considered 
in Ref. [7]. Our results should be contrasted with the 
mean-field analysis of Ref. [7] which favors a gapless SL, 
rather than SPS, in the relevant parameter range. We 
attribute the difference to the fact that the mean-field 
approach [7] neglects essential gauge fluctuations, which 
are accounted for by VMC. 
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FIG. 2: Energies of AFM, ASL, SPS, and dimerized phases 
compared for 10 x 10 system. AFM state is favorable at 
J2/J1 ^ 0.08; ASL and SPS states have energy lower than 
AFM at J2/J1 > 0.08, but their energies are very close. 
Spontaneous breaking of the rotational symmetry occurs at 
J2/J1 0.3, giving rise to the dimerized state. 




FIG. 3: Energies of the SPS state for different values of pairing 
amplitude A as a function of pairing phase 8 compared to the 
ASL energy (for system size 14 x 14). SPS state is favored, 
and its energy is minimized for A ~ 0.1,9 ~ 1.1. 



Spin Hamiltonian We start from an effective spin 
Hamiltonian, the J\ — J2 spin model, 



H = J\ Si ■ Sj 
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where (ij), ((ij)) denote nearest-neighbor and next- 
nearest-neighbor sites. To establish the connection with 
the Hubbard model, we calculate the parameters of the 
spin model from the perturbation theory in (t/U) 2 jS], 
finding, to second order in (t/U) 2 , 
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The SL phase in the Hubbard model then ranges from 
J2/J1 m 0.07 at the antiferromagnetic (AFM) transition 



to J2/J1 ~ 0.12 at the semi-metal transition. Both ex- 
change couplings are antiferromagnetic, thus the effective 
spin model is frustrated. In this study we ignore higher 
order terms in (t/U) 2 as well as third nearest neighbor 
and ring exchange terms [21 . At the relevant U/tm 4.3 
higher order terms are likely too small to affect the results 
and third nearest neighbor terms are non-frustrating and 
will just renormalize the effective J\ . 

Quantum fluctuations are particularly important in 
this model as they are enhanced by the low cooridina- 
tion number of the honeycomb lattice and competition 
between the nearest and next-nearest-neighbor exchange 
interactions. We contrast the phase diagram of the quan- 
tum model (Fig. 1) with that of the classical J\ - J2 
model. The latter model exhibits two phases: Neel anti- 
ferromagnet, with opposite spin polarization on the two 
sublattices of the honeycomb lattice, and (an incommen- 
surate) spiral ordering [13 . The Neel state, which is 
favored by the nearest-neighbor coupling, is realized at 
small J2/J1 < 1/6, where the nearest-neighbor coupling 
dominates. When J2/J1 > 1/6, spiral order with an 
incommensurate wave vector sets in. The origin of the 
spiral order is especially clear at very large !/> 1, when 
nearest-neighbor coupling is negligible, and the two tri- 
angular sublattices can be thought of as independent; the 
ground state on each sublattice is a spiral state in which 
spins on three sub-sublattices are rotated by 2-7r/3 with 
respect to each other. 

We find the quantum fluctuations drastically alter this 
phase diagram. The spiral phase is destroyed at interme- 
diate J2/ Ji, giving way to a disordered spin- liquid phase 
at J2/J1 < 0.3, and a dimerized phase at J2/J1 ^ 0.3. 
The Neel state survives at small J2/J1 0.08, albeit 
with reduced magnetization. 

Ansatze. The main goal of this work is to understand 
the nature of the spin- liquid phase. Toward that end we 
focus on two primary types of wave-functions: (gener- 
alized) Huse-Elser jlOl [T7] states and resonating valence 
bond (RVB) states [5]. 

The former of these is chosen as a good ansatz for the 
AFM state. In the Huse-Elser wave function, the phase 
of the wave function is fixed by the Marshall sign rule 
and for the real part we optimize a separate variational 
two-body parameter C(r) for each unique vector r. 

The RVB state is represented as 
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where {T>} is a (generically non-nearest neighbor) dimer 
covering of the lattice. Different choices for Av corre- 
spond to qualitatively different types of wave functions. 
RVB states are good ansatze for (gapped and gapless) 
spin-liquid states as well as dimer states. 

One approach for selecting these amplitudes is to write 
down a large (but not complete) set of parameters speci- 
fying the RVB amplitudes and then optimize over them. 
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We parameterize Av so as to be able to represent all 
states of the form PYlk( u k + ll /c c l -|- c -fc 4.) |0) for any real 
Uk,Vk where P projects out double occupancy. We call 
these generic RVB states. Optimization for these (and 
the Huse-Elser) states is done via stochastic optimization 
|18j . This parameterization has two drawbacks. First, 
because optimization of a large set of parameters runs 
the risk of being stuck in local minima, we are not guar- 
anteed to find the best state. Secondly, some spin liquid 
states, such as SPS (see below), are not encompassed 
by this parameterization. Because of these concerns, we 
also generate RVB amplitudes in a more physically mo- 
tivated way allowing for fewer parameters that can be 
directly swept over. 

This alternative approach uses the Schwinger fermion 
representation of the spin model combined with 
Gutzwiller projection [6]. In this approach, the 
spin operator on the ith site is related to fermionic 
creation-annihilation operators /, p as follows, Sf = 
^2s s' fists' /*«' ' an d a constraint of one fermion per 
site is imposed, flsfis — L Wave functions of the 
spin model are obtained by Gutzwiller projection of the 
fermionic many-body wave functions, which projects out 
sites with double or no occupancies. 

The fermionic wave functions are then generated as 
ground states of a quadratic Hamiltonian on the honey- 
comb lattice 

a, / x .c./^ -X A -^/;. fir:-; ■'>■<■■ (5) 

{ij),s ij 

which includes nearest-neighbor hopping, and supercon- 
ducting pairing. In principle, a small second neighbor 
hopping should also be included as suggested by the RVB 
mean-field theory of this state but we ignore it in this first 
pass at the problem. The parameters {A} are chosen in 
such a way that the ground state energy of the projected 
wave-function is minimized. 

An important advantage of the Schwinger fermion rep- 
resentation is that there exist simple choices of {A}, with 
just hopping matrix element between nearby neighbors, 
which describe different types of candidate spin liquid 
states that have been proposed. Lee and Lee [T5], and 
later Hermele [15] conjectured the existence of an alge- 
braic spin liquid (ASL) on the honeycomb lattice, which 
is characterized by gapless spin excitations with Dirac- 
like spectrum, similar to that in graphene. This corre- 
sponds to the nearest-neighbor tight-binding model, with 
Aij = for all Very recently, Lu and Ran [7] sug- 
gested two possible gapped spin liquid states, which differ 
by their projective symmetry group [5 ; we will call these 
states s-wave spin liquid (sSL), and sublattice-pairing 
state (SPS). sSL is obtained by considering real A^'s, 
which are rotationally and translationally invariant, for 
sites i, j which are nth nearest neighbors or closer. SPS 
is characterized by complex pairing amplitudes, with op- 
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FIG. 4: Energy difference between SPS state with optimized 
(for 14 x 14) pairing parameters (A ~ 0.1, 8 « 1.1) and ASL 
state as a function of system size. The energy difference ex- 
trapolates to a non-zero value in the thermodynamic limit 
N — > 00. Non-monotonicity of the points is a result of incom- 
mensurability effects with the lattice. 
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FIG. 5: Energy difference for system size 14 x 14 between 
SPS (S = 0.1, 8 — 1.1) and ASL states, as a function of the 
frustration parameter Jij J\. AE > indicate the SPS is 
variationally preferred. The energies of the two phases cross 
at J2/J1 ^ 0.2 and further study is warranted to determine 
if SPS with different parameters still dominates ASL in this 
regime. 
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where i,j are next-nearest neighbors. Within their RVB 
mean-field theory, Lu and Ran find that ASL is favored 
over sSL and SPS when J2/J1 < 0.3 which encompasses 
the entire regime relevant to the spin-liquid in the Hub- 
bard Model. In all the spin-liquid ansatze, the symme- 
tries (translational, time-reversal, rotational symmetry) 
of the honeycomb lattice are respected. 

Results. Using the variational approach, we have 
mapped out the energies of various phases, illustrated 
in Fig. [2] We find a phase transition between AFM and 
SL at J2/J1 ~ 0.08. Both phases are also found in the 
Hubbard model jl], and the transition point is remark- 
ably close to that in the Hubbard model. At higher frus- 
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FIG. 6: Dimer-dimer correlation function, defined as in [4] 
for a SPS state. Green is the reference slice. Red indicates a 
positive correlation with the reference slice and blue a nega- 
tive correlation. The SPS has the same positive and negative 
correlations as the dimer-dimer correlations found in 0j. 

tration parameters (J2/J1 > 0.3), we find the rotational 
symmetry of the RVB states is broken giving a dimerized 
state a lower energy then that of the SL phase. This is 
seen by optimizing RVB amplitudes up to third nearest 
neighbors (with the other amplitudes fixed as in the ASL 
state) and is consistent with findings from exact diago- 
nalization studies on small clusters f I12j ) which have 
suggested dimerized states. Having identified the loca- 
tion for the spin-liquid phase, we turn our attention to 
identifying its nature. 

In establishing the form of spin liquid state, we focus 
on the ASL and the SPS state, which is variationally the 
lowest gapped state we find. We do not consider sSL as 
our optimization over generic RVB states (which includes 
sSL) does not find a lower state then SPS. We notice that 
the energy difference between ASL and SPS is very small, 
and a more careful study is needed to distinguish between 
them. 

To establish whether SPS is more favorable than the 
ASL, we have optimized the SPS energy with respect to 
pairing amplitude and phase. We first consider J2/J1 = 
0.1, and a 14 x 14 system. By mapping out the energy 
as a function of A, 8 (see Fig. pjj, we have established the 
optimal values A/t « 0.1, 8 w 1.1. 

The energy of the SPS state with those parameters 
is lower than that of ASL, suggesting that ASL is un- 
stable with respect to pairing that opens a gap. How- 
ever, the energy gain due to the gap opening is so small, 
that one may doubt whether it survives in the ther- 
modynamic limit. To answer this question, we studied 
finite-size scaling of the energy difference E$ps — Easl 
at the parameter values A/t = 0.1,8 = 1.1 considering 
N = {10, 11, 13, 14, 16, 17, 19}. The result, illustrated in 
Fig. [4] clearly shows that the energy difference extrap- 
olates to a non-zero value in the thermodynamic limit 
1/N 0, indicating that SPS is the ground state. 

We have repeated the comparison between energies for 
ASL and SPS in the whole range of frustration parameter 
0.05 < J2/J1 < 0.25; the result for the 14 x 14 system 



is illustrated in Fig. [5] We found that SPS state is fa- 
vored in the range 0.05 < J2/J1 < 0.2, and at larger 
values of the frustration parameter the energies of the 
two phases are swapped. Additionally at these higher 
frustration parameters, we find a generic RVB state that 
does not break sublattice symmetry and has a lower en- 
ergy then the ASL state. Because we have not studied the 
finite-size effects or optimized carefully the SPS param- 
eters at these higher frustration parameters, this could 
either point to a series of phase transition between the 
different states, or to the SPS gap becoming too small to 
be resolved without more careful optimization and finite 
size extrapolation. Further work is needed to distinguish 
between these different scenarios. 

Having identified SPS as the variationally lowest en- 
ergy state, we examine some of its properties. We start 
by naively estimating the gap. We accomplish this by as- 
suming that the gap corresponds to the BCS pairing gap 
of the tight-binding model obtained by a mean-field de- 
coupling of the Schwinger fermion model. For the 14 x 14 
system we find 



which compares favorably with Hubbard model result [I] , 
AE ss i/13 for the same system size (in the thermody- 
namic limit AE t/40). 

We also look at the dimer-dimer correlation function 
for an AFM and SL state (see Fig. [5) . We find that the 
dimer-dimer correlations of the SL state are positively 
(respectively negatively) correlated on exactly the same 
dimers as the Hubbard model at U/t ps 4.0 |4J. It should 
be noted that although the ASL state has a similar look- 
ing dimer pattern, due to the small gap of the SPS, the 
AFM state looks qualitatively different. 

Discussion. In conclusion, we have studied J\ — J2 
model on the honeycomb lattice, finding evidence for SL 
phase in the range 0.08 < J2/J1 < 0.3. We have accu- 
mulated strong evidence that the SPS state describes the 
spin-liquid phase seen in the Hubbard Model. Beyond 
having the transition happen in almost exactly the cor- 
rect place, we find it to be the variational lowest energy 
state beating out the gapless ASL state. Additionally, in 
naively estimating the gap, we find that it corresponds 
well to the Hubbard gap. Moreover, the dimer-dimer 
correlations closely match those of the Hubbard model. 
We should note that for SU(2), as opposed to large N 
SU(N) spins, it is not clear whether the SPS state repre- 
sents a phase of matter that is distinct from the simplest 
short ranged RVB phase obtained, e.g., in the quantum 
dimer model |22| . 
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